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I. INTRODUCTION 

Since the first experimental observations (TJ El El of Bose- 
Einstein condensation in dilute alkali gases, much interest has 
been directed to the study of quantum vortices in such sys- 
tems, both experimentally El and theoretically (for re- 
views see ||7l|8l[9l[lQl|). On one hand, dilute condensates of- 
fer an opportunity to gain insight into the physics of superflu- 
ids and vortices familiar from condensed-matter systems. The 
weakly interacting nature of dilute condensates makes them 
amenable to tractable theoretical approaches, complement- 
ing the precision of control and observation aflTorded by dilute 
atomic-gas experiments |5|. However, there are also marked 
diflferences between the vortex physics of dilute condensates 
and that of bulk superfluids. In contrast to bulk superfluids, di- 
lute condensates have inhomogeneous density profiles, lead- 
ing to new vortex physics 01111 US- Furthermore the behav- 
ior of dilute gas condensates is in many cases subtly diflferent 
to that expected from experience with bulk superfluids: exam- 
ples are the much slower collision rates in the dilute conden- 
sates, resulting in the breakdown of two-fluid hydrodynamics 
|[T3l[T4]|, and the essential role of dynamical instabilities in the 
formation of vortex lattices in stirred condensates 1151 . 

Early numerical investigations by Butts and Rokhsar |[T6]| 
revealed some of the the nature of the rotating condensates 
at zero temperature. As in bulk superfluids 11711 , dilute con- 
densates seek to mimic rigid-body rotation ifTSl by acquiring 
vortices in a regular array. Condensates adjust their rotation 
both by undergoing discontinuous phase transitions (161 [191 
between rotational-symmetry classes by nucleating additional 
vortices, and also by adjusting their vortex distributions within 
a symmetry class. At a critical imposed rotation frequency, 
the ground state of the system makes a transition from an ir- 
rotational, vortex-free state to a central- vortex state |20|, with 
angular momentum per atom i = {Lz)INfi = 1. This tran- 
sition is continuous, and the system possesses a continuum 
of oflT-axis vortex states with < ^ < 1 . Unlike the stationary 
^ = and £ = I states, these intermediate states are not ground 
states in any rotating frame lfT6l . and are thermodynamically 
unstable in all frames. However, they are dynamically sta- 



ble, and at zero temperature the vortex follows a stable orbit 
about the trap center. Such states represent the simplest non- 
trivial example of vortex dynamics in BEC. The precession of 
such a vortex can be viewed as the vortex being carried by its 
own self-induced superflow 1211 . arising due to the eflTective 
boundary conditions imposed on the superflow by the inho- 
mogeneity of the condensate orbital density |22|. Corrections 
to this picture due to the finite extent of the vortex core have 
been investigated in [23 1 and |24|. 

Further complexity is attained in finite-temperature scenar- 
ios, where the vortex may be subject to additional forces re- 
sulting from its interaction with the thermal field 1 25 , 26] . If 
the angular velocity of the thermal cloud is diflTerent from that 
of the vortex, the vortex experiences a frictional force which 
causes it to drift radially. For thermal-cloud angular velocities 
below some threshold |7|, this interaction induces the expul- 
sion of the vortex from the condensate (271 . The dynamics 
of a vortex at finite temperature are of interest both as a test 
of theories of finite-temperature condensates ll28l[29ll3Qll3n , 
and for understanding the nature of vortex interactions with 
bulk thermal flows, which has proved a controversial issue 
1251 l26l [32l . In this paper we realize finite-temperature pre- 
cessing single-vortex states as equilibrium configurations of 
a conserving (Hamiltonian) classical-field system |31|, with 
fixed internal energy and angular momentum. Such Hamil- 
tonian classical-field methods and closely related stochastic- 
field methods | 33 , 34] have proven useful both in equilibrium 
and nonequilibrium scenarios (IT! l35l [36l ITTl [38l [39l l4Ql. 
Working at fixed angular momentum allows us to create 
rotational-symmetry-broken precessing- vortex states at finite- 
temperature equilibrium of the classical field. In such configu- 
rations, the vortex is at rotational equilibrium with the thermal 
cloud, and the dissipative force responsible for vortex decay 
lETlEnmil vanishes. We extract the symmetry-broken con- 
densate orbital from the classical-field equilibrium, and ob- 
serve both the Goldstone mode associated with breaking ro- 
tational symmetry, and the filling of the vortex core by the 
thermal component of the field. Increasing the field energy 
at fixed angular momentum we observe an increase in vortex 
precession radius and frequency, ultimately leading to the ex- 
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pulsion of the vortex from the condensate and the decoupHng 
of the condensate rotation from that of the thermal cloud. 

This paper is organized as follows: In Sec.|Il|we introduce 
the formalism we use in this work, and discuss its interpreta- 



tion. In Sec. Ill we discuss our simulation procedure, intro- 
ducing in turn the parameters of our system, and the details 



of our numerical approach. In Sec.|IV]we present results of 
representative simulations, and discuss techniques used in the 
extraction of physical information from the simulation trajec- 
tories. In Sec. |V] we discuss the dependence of the system 



behavior on the internal energy of the field, and in Sec. [VI 
summarize and present our conclusions. 



we 



II. FORMALISM 
A. Review of Hamiltonian classical field method 

In this paper we employ classical-field-theory techniques, 
which have recently been reviewed at length in [31 1. Here we 
recapitulate the basic elements of the Hamiltonian classical 
field method, and briefly discuss the appropriate interpretation 
of its application here. 

At its simplest level, classical-field theory results from the 
replacement of a second-quantized (Bose) field Hamiltonian 
with its classical analogue. This is motivated by the observa- 
tion that in regimes of thermal behavior, the dominant charac- 
teristic of the atomic Bose-field system is its multimode, self- 
interacting nature. Indeed, under conditions of high mode oc- 
cupation, field commutators become relatively unimportant, 
and a satisfactory description of the bosonic field can be ob- 
tained from a classical-field model. More generally, the dy- 
namics of the classical-wave system itself are both of intrin- 
sic interest B2l . and able also to provide qualitative insight 
into the dynamics of degenerate Bose-gas systems in scenarios 
where a precise identification of the method with the many- 
body field theory is impractical |38|. In the present work we 
consider a closed, Hamiltonian system, with fixed angular mo- 
mentum. This allows us to characterize the motion of a quan- 
tum vortex at equilibrium with a rotating thermal component 
of the field. 

The essential ingredient of the classical-field method as de- 
veloped and refined by Davis, Blakie and co-workers lISTllS?! 
[35l|43 | is the projection operator. While the classical dynam- 
ics of infinite-dimensional systems is able to be perfectly well- 
defined (see e.g. 1441 ). the projection operator restricts the 
system to a finite number of degrees of freedom [45 1 . This 
projection reduces the dynamics of the field to a finite Hamil- 
tonian system, suitable for numerical implementation. Fur- 
thermore it acts to regularize the divergences resulting from 
the contact scattering potential assumed in the field theory 
1461 [47 1, such that a formal renormalization can then be per- 
formed if so desired BSl . 

The finite-dimensional Hamiltonian system so obtained ex- 
hibits ergodic behavior | 49 1, and this provides a useful tool for 
analyzing the equilibrium properties of Bose-gas systems. By 
time-averaging functionals of the projected field along a suffi- 
ciently long field trajectory, one expects to obtain their mean 



values in the microcanonical ensemble of field configurations 
satisfying the constraints of the conserved first integrals of the 
system |50|. These mean values then serve as estimates of the 
analogous quantum- statistical-mechanical correlations of the 
many-body quantum system. With a careful choice of pro- 
jector, the classical (equipartitioned) equilibrium distribution 
can be chosen to be coincident with the high-occupation limit 
of the true bosonic system, allowing for highly accurate esti- 
mates to be made for the quantum system |[36ll . 

The utility of the classical-field method in its application 
to dynamical systems is its ability to include the eff'ects of 
such thermal behavior nonperturbatively in the description of 
the 'macroscopic' dynamics of the field |[38l[39l. The coher- 
ent fraction's evolution is influenced by the (locally) thermal 
behavior of the incoherent fraction, including both mean- field 
repulsion and dissipative eff'ects in its resulting dynamics. The 
coupled dynamics of the two components are described also 
in the method of 1 13|, with the crucial diff'erence that in the 
classical-field approach no distinction between condensed and 
thermal material is made in deriving equations of motion or 
implementing their solution numerically. Instead the conden- 
sate is identified a posteriori from the fluctuation statistics of 
the field. 

In this work, while we wish to study the behavior of the 
classical-field system at equilibrium, there is an essential dy- 
namical aspect to the equilibrium, i.e., the precession of the 
vortex. Our interest is in the dynamics of the vortex, and thus 
of the condensate, in the presence of the thermal component of 
the field. We thus follow 1 39 1 in employing the classical-field 
method to include the eff'ects of the thermal component of the 
field on the motion of a collective excitation of the conden- 
sate. In contrast to that work the collective excitation we study 
here (the precession of the vortex) is stable at thermal equilib- 
rium of the system. The equilibrium trajectory of the vortex, 
which is determined by the balance of coherent and thermal 
forces acting on it, defines the frame in which the rotational- 
symmetry-broken condensate mode is stationary. The ergodic 
character of the field evolution tends to restore this broken ro- 
tational symmetry as time progresses, and we are thus lead to 
abandon formal ergodic averages and consider the coherence 
properties of the field on short time scales, as we discuss in 
SecHVl 



B. System 

We consider here a system harmonically trapped with a 
trapping potential isotropic in the xy plane: 



y(x)= ^[a;?(x2+/) + a;,V]. 



(1) 



We choose highly oblate harmonic trapping, with trapping fre- 
quency sufficiently high that no modes are excited in the 
z-direction. A two-dimensional (2D) description is thus valid 
provided that 1,5 1.L 



(2) 
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where ji and T are the system chemical potential and temper- 
ature respectively (see also the discussion in Ref. |38 1). Re- 
stricting the system to 2D permits the dynamics of a point 
vortex interacting with the thermal component of the field to 
be studied, while removing complexities such as line bend- 
ing and Kelvin-wave excitations of a vortex filament of finite 
extent, which provide additional mechanisms for vortex dissi- 
pation |52|. Following Ref. |38|, we choose physical parame- 
ters corresponding to ^^Na atoms confined in a strongly oblate 
trap, with trapping frequencies {ojr.ojz) = 2;: x (10,2000) 
rad/s. The s-wave scattering length is a = 2.75nm, placing 
our system in the quasi-2D regime (/^ = ^Jf^JnuOz » a), with 
an eff'ective 2D interaction parameter [/2D = 2 ^|27lf^a|mlz. In 
our simulations we will often compare our systems to what we 
will refer to as our principal ground state, which we take as the 
ground Gross-Pitaevskii (OF) eigenmode | 20 | of the trapped 
system in a frame rotating at angular velocity Qq = 0.356l)^, 
with eigenvalue (^ig)^^ = lOHoJr in that frame. This state 
consists of A^o = 1.072 x 10^ atoms, and contains a sin- 
gle on-axis vortex carrying unit angular momentum per atom 
(angular momentum and rotation will always refer to that 
about the z axis in this paper). In an inertial (non-rotating) 
frame this state has eigenvalue yUg = 10 .SSfioJr, and energy 
Eg = 7.646 X lO'^hoJr. In this paper we will refer generically 
to any such inertial frame as a laboratory frame (lab frame). 



C. Equations of motion 



To derive the classical-field equation of motion, we begin 
with the second-quantized Hamiltonian of a trapped Bose gas 
in the ^-wave regime 

H = Jd^x ^\x)H,p^(x) -^^ Jd^x ^\x)^\x)^(x)^(x), 

where ^(x) is the bosonic field operator satisfying 
[^(x), 4'^(x')] = ^(x-x'), and the interaction is assumed to be 
a contact potential with U = Anfi^alm, a the scattering length 
and m the atomic mass. We choose to write the Hamiltonian in 
a frame rotating at angular velocity about the z axis. As the 
interaction potential is rotationally invariant this choice aff'ects 
only the definition of the single-particle Hamiltonian 



(4) 



We introduce a cutoff' energy Er, delimiting a low-energy re- 
gion (L) of the system consisting of single-particle modes 
[eigenmodes of Eq. Q] with single-particle energies 6„ < 
Er, which is sometimes referred to as the condensate hand 
1311 . The remaining modes comprise the complementary 
high-energy region (noncondensate band H). The cutoff' Er 
is chosen to be at a sufficiently high energy that the interact- 
ing Hamiltonian Eq. ^ is approximately diagonal at the cut- 
off* 13 3 J . This division is enforced formally by the projection 



operator defined 

^/(x)^^0,(x) r JV0:(y)/(y), 



(5) 



where summation is over all single-particle modes satisfying 
€n < Er. We refer to the number of such modes spanning the 
low-energy region as the condensate-band multiplicity (At). 
We use this projector to define the projected field operator 



iJf(x) = Y,^nM^) = P^(x). 



(6) 
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We choose the cutoff' such that Er < fioj^, and thus all ex- 
cited z-axis modes are excluded from the low-energy basis 
over which i//(x) is expanded. A classical theory is obtained 
by demoting the operators an in Eq. ^ to classical vari- 
ables an, thereby defining a projected classical field i//(x) = 
ZneL (^n(pn(^)- The Corresponding classical-field Hamiltonian 
is 



HcF = J d^x 



r(x)ifsp^(x)+^|iA(x)|4, (7) 



which is formally a classical Hamiltonian with 2M canon- 
ical degrees of freedom Qn} related to the field vari- 
ables {Qf„,Qf*} by a canonical transformation 1531 . It then 
follows that the Hamilton's equations can be expressed in 
terms of projected functional diff'erentiation of the classical- 
field Hamiltonian 1331 as 



dilfix) SHcF 



dt 



Si//*(x) 



p[{H,^^U2Dmf)m}, (8) 



which is the projected Gross-Pitaevskii equation (PGPE) (see, 
G-g-. 11311 ). The conservation of field energy E[i//] = Hqf and 
normalization N = J d^x\i//(x)f' under the action of the PGPE 
follow immediately from the nature of ^cf as a time-invariant 
classical Hamiltonian, and its invariance under the transfor- 
mation </^(x) il/(x)e^^ respectively (441 . Furthermore the 
rotational invariance of the Hamiltonian ^cf implies that the 
angular momentum 



r(x)L,i/f(x) 



(9) 



[where the bar denotes spatial averaging: A = 
J d^xi//* (x)Ai//(x)] is also conserved by the evolution of 
the PGPE. 



III. SIMULATION PROCEDURE 
A. Numerical implementation 

1. Dimensionless units 

For the purpose of the implementation of the PGPE, it is 
convenient to express the physical quantities {r, oj, E, t} in a 
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dimensionless form, which we write as {f,a),E,t}. These 
quantities are related by the expressions r = fro, oj = cocor, 
E = Ehojr, t = toj~^, where the radial oscillator length 
ro = ^JfilmcOr- In reporting our results we will often refer 
to times in units of the trap cycle (abbreviated eye), i.e., the 
period of the radial oscillator potential Tosc = InjcOr- 



2. PGPE implementation 



Here we briefly discuss our implementation of the PGPE 
in terms of rotating-frame harmonic-oscillator states, as de- 
scribed in |[34j|38l. In dimensionless form the PGPE [Eq. ([8])] 
becomes 

= ^{[ - y + y + i^de + ^|^Al']^|, (10) 

where the dimensionless efl'ective interaction strength A = 
U2'Dlrpi(ji)r- We proceed as in |[34l by expanding ij/ as 



(11) 



{«,/} 



over the Laguerre- Gaussian modes 



(7r(^ + |/|)! 



^//^p|/|^-PV2^|/|(p2)^ (12) 



which diagonalize the single-particle (non-interacting) Hamil- 
tonian H[ 



sp 



V2 



H%Yni{h ^) = [- y + -+i^de]Yni{f, 6) = E^iYniif, e\ (13) 



with eigenvalues E^^ = 2^ -h |/| - Q/ -h 1. The energy cutofl' 
defined by the projector V requires the expansion of Eq. ( pT) 
to exclude all terms except those in oscillator modes ?„/(r, 6) 
for which E^^ < Er, i.e. those modes satisfying 



2n-\-\l\-ai-\-l < Er. 
This yields the equation of motion 

= E^fnl + AFnl(i//). 

where the projection of the GP nonlinearity 

^2n 



(14) 



(15) 



r*Z7T r^oo 

Fnm= do rdrY:,(r,emr,efi/f(r,ei (16) 
Jo Jo 

is evaluated by the technique presented in |[34l . 



B. Microcanonical evolution 



1. Microcanonical formalism 



Davis, Blakie and co-workers have shown |[35j|54l that the 
PGPE can be used to evolve a general configuration of the 
classical field to a thermal equilibrium. In contrast to those 
works, in which the initial field configurations were parame- 
terized solely by their energies 



EW = HcfW 



(x)//sp(A(x)+-^|<A(x)|\ (17) 



in this work we also fix the conserved angular momentum 



L[^] 



(18) 



to a prescribed nonzero value |55|. We therefore expect that 
the field evolves under the action of the PGPE to a ther- 
mal equilibrium consistent with the choices of conserved en- 
ergy [Eq. ([17])] and angular-momentum [Eq. ([18])] first inte- 
grals. In this paper we choose L[i//] = N[i//]fl = NqH = Lq, 
and thus explore the manifold of classical-field configurations 
with unit angular momentum per particle. At zero temperature 
(corresponding to complete condensation in the classical-field 
model employed here, in which no representation of quantum 
depletion is included), the corresponding microstate (unique 
up to a choice of phase) is our principal ground state, con- 
taining a singly charged, on-axis vortex, with energy E[i//] = 
Eg. By choosing initial configurations with E[i//] > Eg and 
Lzii/^] = Lo, we investigate the thermal equilibria that result 
when atoms are thermally excited out of the ground conden- 
sate mode, while conserving the total angular momentum. 



2. Choice of frame 

In applying our classical-field method, we must make a 
choice of the rotation rate of the frame in which the PGPE is 
defined, which we will denote by Q.^ in the remainder of this 
paper. The only non-trivial consequence of this choice 15611 
is the precise definition of the projector P = P{Er\ Qp), and 
thus of the set of modes which span the low-energy space L. 
This dependence of the projector on Qp generalizes the famil- 
iar concept of an energy cutoff* in classical-field theory. As is 
well known [31 ], the thermodynamic parameters of the clas- 
sical field at equilibrium depend both on the conserved first 
integrals of the system and the imposed cutoff*. In the present 
work, we expect that the equilibrium of the system minimizes 
the free energy 



F = E -TS -aL, 



(19) 



where the thermodynamic angular velocity = (dE/dL)s 
(571. The utility of the rotationally invariant (Gauss-Laguerre) 
projector used here is that the classical field evolution is not 
dynamically biased towards a particular rotating frame [38] . 
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i.e., we retain a (generalized) Ehrenfest relation iSE 



dLy 



(20) 



for the angular momentum L^. This ensures that for the 
isotropic trapping we consider here, angular momentum is 
conserved, and further that the field is free to rotate at the 
angular velocity D which minimizes Eq. ([19]), which will in 
general not be equal to the projector angular velocity Qp. The 
relationship between the first integrals {E, L) and the thermo- 
dynamic variables (r, Q) of course depends on Qp, as diff'er- 
ent choices of cutoff' yield diff'erent Hamiltonian systems. We 
demonstrate this dependence by explicit calculation in Ap- 
pendix [A| and note that the use of a harmonic basis allows the 
PGPE theory to be extended to include a mean-field descrip- 
tion of above-cutoff' atoms, whereby insensitivity to (moder- 
ate) changes in the cutoff' is regained (see |31 1). However, as 
a first application to the precessing- vortex system, we restrict 
our attention to a strictly Hamiltonian "PGPE system" l35l , 
with low-energy region L defined by application of the pro- 
jector in the frame in which the ground state was formed (i.e., 
we take Ilp = Qq = O.SSci;^), and investigate the behavior 
of this particular Hamiltonian system by varying its energy at 
fixed angular momentum. 



3. Procedure 

For almost any initial configuration far from equilibrium, 
we expect the field to approach equilibrium over time B9l , 
due to the ergodicity of the system. We thus form our initial 
states as randomized field configurations under the projector 
defined by cutoff' energy Er = SQigJa^ = SOHaJr and rotation 
frequency Qp = 035<jJr, with the same values of normaliza- 
tion and angular momentum L as our principal ground state, 
but with higher energies. These states are essentially formed 
by 'mixing' a randomized high-energy state with the princi- 
pal ground state in a procedure similar to that described in 
II3TII . The initial-state information will persist in the precise 
details of the ffuctuations | 42 1 since the field equations are de- 
terministic, although reversibility of the evolution is lost after 
some time to numerical error |49|. The values of the con- 
served first integrals L and E are fixed to the desired values to 
within specified tolerances (AL/Lq = lO'^A^/^g = lO'^). 
We form initial states in this manner with a range of energies. 
In each case the central vortex of the principal ground state 
remains, however the density profile of each state is severely 
distorted. As all our states have L = NqH, the inertial-frame 
and projector-frame energies are in each case related by a con- 
stant shift resulting from the fixed rotational kinetic energy of 
the field: (E)o = (E)q^ -h fi^^No. This shift is of order 5% for 
the ground-state wavefunction. In this paper we will discuss 
the energies of states in the lab frame, and quote them as a 
multiple of the ground- state energy Eg in that frame. 

Having formed the initial states we evolve them in the pro- 
jector frame for a period of 10"^ trap cycles, with an adap- 
tive integrator [59J with accuracy chosen such that the rela- 



tive change in field normalization is < 4 x 10"^ per time step 
taken |38|. By 9 x 10^ trap cycles the states have reached 
equilibrium, as evidenced by the settling of the distribution of 
particles in momentum space |[54l [601, and we perform our 
analysis of steady-state properties on the final 10^ trap cycles 
of evolution data. 



IV. RESULTS 

After the equilibration period, the behavior of the classi- 
cal field simulations is qualitatively similar for energies in the 
range £" = [1.04, 1.15] £g. Each exhibits a central high-density 
region containing a vortex which is displaced off' and precess- 
ing about the trap axis. Viewed in the lab frame, the preces- 
sion is in the same sense as the vortex (i.e. positive rotation 
sense, corresponding to angular-momentum numbers m > 0), 
and is stochastic in nature, with the vortex position ffuctuating 
about a circular mean orbit, in the presence of strong density 



ffuctuations of the background ffeld. As noted in Sec. IIIB 2 



the frequency of the vortex precession is in general not equal 
to the angular velocity of the projector defining the micro- 
canonical system. The condensate is located entirely in the 



central bulk (see Sec. IVB ); outside this region the field ex- 
hibits the high-energy turbulence 1 38 1 characteristic of purely 
thermal atoms in the classical-field model. In Fig. [TJa-f) we 
plot classical-field densities spanning a single orbit of the vor- 
tex in a simulation with conserved energy E = 1.10£'g. The 
behavior presented in Fig. [T] is representative of all simula- 
tions in which the equilibrium condensate contains a precess- 
ing vortex. From a range of simulations we observe that the 
radius at which the vortex precesses increases as the energy 
of the field is increased. At higher energies it approaches the 
edge of the condensate and dXE ^ 1.16£^g it is ultimately lost 
into the violently evolving peripheral region of surface exci- 
tations and short-lived phase defects. As the field energy is 
decreased, the vortex approaches the trap center. However, 
we expect that as the energy of the field is lowered towards 
the ground state, the excitation becomes too weak to reliably 
generate thermal statistics on the time scales of interest, and 
we suspect that for low energies the background field 'seen' 
by the vortex may not appear thermal on the time scale of 
its precession. Indeed in such a low temperature regime the 
eff'ects of quantum ffuctuations neglected in our model will 
become important. We therefore exclude simulations with en- 
ergies E < 1.04£^g from our analysis. A detailed discussion of 
the dependence of system observables on the internal energy 
of the ffeld is presented in Sec.[V| 



A. Vortex precession 

In order to characterize the vortex motion, we track the 
vortex location as observed in the laboratory frame. In prac- 
tice we recorded vortex locations in the simulation with E = 
1.10£^g, at a frequency of 25 samples per cycle for a period 
of 400 cycles, starting at ^ = 9000 eye. We ffnd that the dis- 
placement of the vortex from the trap center ffuctuates about a 
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mean value of r\ = 2.04ro, and that the series of displace- 
ments measured (from the trap center) has standard devia- 
tion (Tr^ = 0.09ro, which is of same order as the extent of 
the vortex core (which is approximately given by the healing 
length T] = 0.20ro |20, 61 1). The normalized distribution of 
displacement radii measured over the 400 eye. period is dis- 
played in Fig.[2ja). The plotted distribution shows some neg- 
ative skew but this is not a consistent feature across the sim- 
ulations performed. To determine temporal characteristics of 
the vortex trajectory, we define a complex vortex coordinate 
Zy(t) = Xy(t) -\- iyy(t), where the coordinate pair (xy(tj),yy(tj)) 
specifies the vortex position at time tj, and we calculate the 
power spectrum of this time series. As direct periodograms 
are subject to large variances in the estimates they provide 
at each frequency (62] we split the 10"^ -sample time series 
into 10 consecutive segments of equal length and average the 
power spectra obtained from these shorter series (Bartlett's 
method f63l). In Fig. |2jb) we plot the power spectrum of 
the vortex coordinate Zy over the (angular) frequency range 
oj e [-12.5, 12.5]ajr permitted by the sample frequency of 
our data ||64l. The most conspicuous feature of the spectrum 
is the sharp peak [note the logarithmic scale of Fig.[2jb)] cor- 
responding to the precession of the vortex about the trap axis. 
The next most pronounced features of the power spectrum are 
two clearly resolved spikes at frequencies oj = ±lcL>r- We 
identify oj = lojr sls the universal frequency of the dipole- 
oscillation mode (Kohn mode) of the harmonically trapped 
gas 1 65 1, which is immune to thermalization 1661 . The fre- 
quencies ±loJr thus result from a small dipole oscillation of 
the classical field as a whole. Smaller features can be seen in 



the wings of the power- spectrum distribution, which we asso- 
ciate with thermal fluctuations of the classical field (c.f. dis- 
cussion in (381). From the polar form Zyitj) = |Zv(0)k^^'^^^^^ of 
the vortex coordinate, we observe that the vortex phase 6y does 
not increase linearly with time, but fluctuates randomly as the 
vortex moves in the dynamically evolving thermal field. Af- 
ter some time the vortex phase diff'uses, so that the rotational 
symmetry of the system is restored in the ergodic density of 
classical-field configurations. This is illustrated in Fig. |2jc), 
where we plot the distribution of vortex phases 6y as measured 
in a frame rotating at frequency = 0.26 lOci;^. This fre- 
quency corresponds to the peak presented in Fig.|2jb), and has 
been chosen because it yields the minimum variation in vor- 
tex phase over the 400 eye. period we consider. Histograms of 
the normalized vortex-phase distribution over the first 50 cy- 
cles of this period [blue (black) bars] and over the full 400 eye. 
period [red (gray) bars] are plotted. We observe that at short 
times the fluctuations in Oy are small and reasonably Gaussian 
in nature. Over longer time periods the vortex phase drifts 
significantly. This diff'usion of the vortex phase produces the 
distinctive Lorentz-like shape of the power- spectrum distribu- 
tion plotted in Fig.[2jb). In Sec. IV B 2 we discuss the efl'ects 
of this difl'usion on the 'one-body' correlations of the classical 
field obtained from time-averaging the classical-field trajec- 
tory, and the implications for the identification and character- 
ization of the condensate in these simulations. 
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FIG. 2: (Color online) Quantities characterizing the motion of the 
vortex, (a) Normalized histogram of measured vortex displacement 
magnitudes. The vertical dashed line indicates the mean, (b) Power 
spectrum of the complex vortex coordinate Zy = + iy^, revealing 
the precession of the vortex as viewed in the laboratory frame, (c) 
Normalized histograms of vortex phases 6^ measured in a rotating 
frame (see main text) over periods of length 50 eye. and 400 eye, 
starting from t=9000 eye. 



B. Field-co variance analysis 

We turn now to the issue of identifying the condensate in the 
classical-field solutions. In vortex-free 1671 equilibrium sce- 
narios a method of condensate definition based on sampling 
of the classical-field correlations | 35 1 has proven successful. 
The method, which is in obvious analogy to the Penrose- 
Onsager (PO) definition |14] [681 of Bose-Einstein condensa- 
tion in terms of the (quantum) one-body density matrix, in- 
volves constructing the covariance matrix 



Pij = {c*Ci}t 



(21) 



from the propagation-basis coefficients (Sec. Ill A 2[ ) of a sin- 
gle classical-field trajectory, where the the indices i, j each 
index a quantum number pair (n, /), and denotes an av- 
erage over time. Ergodicity of the field evolution [ 49 1 implies 
that averages taken over sufficiently long times approach av- 
erages over the microcanonical ensemble of field configura- 
tions, in which the coherence of a classical-field condensate 
\4T\ manifests as a single dominant eigenmode of pij. How- 
ever, this method of condensate identification has some sig- 
nificant ambiguities in more general scenarios, in particular 
when symmetries of the Hamiltonian are broken by the ap- 
propriate GP state |[69l . For example, the standard PO crite- 
rion itself may prove inappropriate in situations in which the 
condensate exhibits center-of-mass motion fTOl, prompting 
workers to introduce revised definitions of the density matrix 
appropriate to to such systems, in which the center-of-mass 
motion is eliminated fTTl, l72l . 

Issues of condensate definition also arise in the presence of 
vortices. In scenarios in which several vortices exist, breaking 
the rotational symmetry of an otherwise rotationally invariant 
many-body Hamiltonian, an uncountably infinite degeneracy 
of the one-body density matrix appears, even in the idealized 
Gross-Pitaevskii limit of a zero-temperature system L73J . In 
this paper we have chosen one of the simplest scenarios of 
rotational symmetry breaking in a finite-temperature classical 
field: the precession of a single vortex. Because the presence 
of the vortex distinguishes diff'erent rotational orientations of 
the condensate, we find that in order to determine the conden- 
sate mode we must construct time averages in an appropriate 
rotating frame. The fact that the phase of the vortex diffuses 
over time in any such uniformly rotating frame means that 
we must abandon the notion of ergodic reconstruction of the 
microcanonical density by time averages. However, we ex- 
pect intuitively [35 1 that in order to quantify the condensa- 
tion in the field, the averaging time need only be long enough 
to distinguish the Gaussian ffuctuations of thermal or chaotic 
modes |[74l [tH from the non-Gaussian statistics of a quasi- 
uniform phase evolution characteristic of a condensed com- 
ponent. We find therefore that we can successfully quantify 
the condensation of the field from such short- time averages, 
while exploiting the ergodic character of the evolution over 
longer time scales. In the remainder of this paper, we will use 
the term Penrose-Onsager (PO) procedure generically to refer 
to the construction of the covariance matrix (density matrix) 
Eq. ( [2T] ) by a time- averaging procedure, the precise details of 
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which will be in each case specified in the context. 



1. Lab-frame analysis 

We consider a simulation with energy E = 1.05£^g, which 
exhibits a vortex precessing at a frequency ojy ^ 0.23ajr, at 
a radius r\ ^ 1.3ro, which should be contrasted with the ex- 
tent of the central density bulk rb ~ 5ro. The vortex in this 
case is displaced from its axis by significantly more than the 
extent of its core (healing length r] = 0.20ro), while the core 
remains comparatively close to the center of the region of sig- 
nificant density. We proceed by forming the covariance matrix 
by averaging over 2501 equally spaced samples of the lab- 
frame representation of the classical-field evolution, between 
t = 9000 and 9100 eye. We diagonalize the density matrix to 
obtain its eigenvectors, denoted by Xi(^) (the eigenmodes of 
the one-body density operator) and their eigenvalues, denoted 
by Hi (which are the mean occupations of the modes during 
the sampled period). The index / ranges from to Al - 1 in 
order of decreasing occupation. In Fig. [Sj^a-e) we plot the 5 
most highly occupied eigenmodes of pij, which together con- 
tain some 95.4% of the total field population. The density and 
phase of the most highly occupied mode |;^o(x)] is displayed 
in Fig. [3ja) and contains a singly positively charged central 
vortex (phase circulation of In around its azimuth direction). 
The second most highly occupied mode |;^i(x). Fig. |3jb)] is 
a rotationally symmetric mode with uniform phase, which 
forms a pronounced peak at the origin. The third most highly 
occupied mode Ixii^), Fig.|3jc)] is another vortex mode, with 
vanishing central density. However, this mode is doubly (pos- 
itively) charged, with the phase varying by An around its az- 
imuth. 

We note that these three modes bear a strong resemblance 
to a single vortex GP-eigenmode 0o and the u and v* spinor 
components |76| of its so-called anomalous Bogoliubov ex- 
citation (also referred to as the lowest core-localized state by 
some authors | 77 1). We recall that this is the lowest-energy ex- 
citation of the single- vortex state in the Bogoliubov approxi- 
mation, and its negative energy in the laboratory frame signals 
the thermodynamic instability of the vortex state in that frame. 
The angular-momentum numbers ruu^v = (0, 2) of the particle 
(u) and hole (v) functions result naturally from the transpo- 
sition of the Bogoliubov pairing ansatz to the m = 1 vortex 
scenario |78, 79 1. The collective excitation of the anomalous 
mode results in the displacement of the vortex off'-axis in a di- 
rection determined by the phase of the Bogoliubov excitation 
relative to the condensate (see LSOJ ) in the linear expansion 



i// = ^000 + bu-\- b*v*. 



(22) 



In the (zero- temperature) Bogoliubov approximation the en- 
ergy of the anomalous mode determines the frequency of pre- 
cession of a displaced vortex under its own induced velocity 
field, in the appropriate linear-response limit |7|. However, 
this result has been shown not to hold in the finite-temperature 
self-consistent mean-field theories which generalize the Bo- 



goliubov description to finite temperatures (29]. 

Returning to our classical field simulation, we calculate 
the time-dependent expansion coefficients of the eigenmodes 
Xi(^) over the 100 eye. analysis period, defined as 

atit)^ J dxx;(xMx,t). (23) 

The average values <|a/p)^ of course yield simply the mode 
occupations which are plotted for the 10 most highly oc- 
cupied modes in Fig.|4ja). From the coefficients {a J we form 
the classical correlation functions 13511821 



(2) (i^^.n. 



,2' 



(24) 



which indicate the coherence of the modes over the analysis 
period. We recall that the (local) correlation functions ^^"^ 
adopt values of g^"^ = nl and g^"^^ = 1 for purely chaotic or 
thermal fields and purely coherent fields, respectively 1135117411 . 
The distribution of g^^^ values over occupation numbers is 
shown in Fig.[3];f). We find that the g^^^ values calculated for 
the modes ;^/(x) : </< 4 are all < 1.1. By contrast the 
vast majority of the remaining modes (/ > 4) have g^^^ ^ 2. 
This suggests that these 5 most highly occupied modes are 
coherent, while the majority of modes are to a large degree 
incoherent with respect to the averaging procedure employed 
here. In |38|, the temporal power spectrum of the classical 
field was used to analyze the coherence in the classical field. 
Here we apply a similar technique to extract information from 
the mode coefficients aiif). In Fig.[4jb) we plot the spectra ob- 
tained from these samples, each of which was formed by av- 
eraging spectra estimates from 5 consecutive time series (see 



Sec. IV A ). We observe that the spectrum of each mode coeffi- 
cient presented in Fig.Qb) exhibits several peaks, the highest 
of which is in each case several orders of magnitude greater 
than the others. We interpret the largest peak in each case as 
an indication of the quasi-uniform phase rotation associated 
with coherence in the classical field. We also note that the 
frequency peaks of modes Sindxi are spaced evenly about 
the peak of xo- Specifically we find for the frequencies of 
peak power (ojp^o,c0p,i,ojp,2) = (10.35, 10.10, 10.60)6i;^, corre- 
sponding to a quasiparticle energy of Sq = h(<jjp^i - ojp^o) ^ 
-0.25hojr, in agreement with the vortex precession frequency 
ojy ^ 0.23ojr (determined by the technique of Sec. |IV A[ ) to 
within the frequency resolution Aoj = 0.056i;^. 

A precessing- vortex GP state can of course be understood 
as a linear combination of such angular-momentum eigen- 
modes. In a frame rotating at the vortex precession frequency 
(Qy), the phases of the various components all rotate at a 
single common frequency (the condensate eigenvalue ji/h). 
In the laboratory frame these components have frequencies 
ojm = j^/h -\- flym, and their consequent dephasing over time 
leads to the apparent fragmentation observed in the lab frame 
ll83l . Simulations of classical fields with higher energies, 
which contain vortices precessing at larger radii, when sub- 
jected to the above analysis, again yield a decomposition into 
angular-momentum eigenmodes. However, in such cases we 
find that the most highly occupied mode has m = 0, as the 



9 



(a) I (b) _ ; (c) 




Occupation Ui 



FIG. 3: (Color online) (a-e) Density and phase of the five most highly occupied eigenmodes of the covariance matrix, for case E = 1.05^^ 
(f) Second-order coherence functions of the covariance-matrix eigenmodes versus mode occupation. 



anomalous particle component u has grown, pushing the vor- 
tex further from the trap center and replacing the m = 1 vortex 
mode as most highly occupied. 



2. Rotating frames 

The above results of a naive application of the PO proce- 
dure in the laboratory frame suggest that the classical field 
condenses in a rotating frame. The mean phase rotation of 
the various angular-momentum components is such that upon 
transformation to an appropriate rotating frame these modes 
collectively exhibit quasi-uniform phase rotation at a common 
frequency. We therefore expect that transforming the field to 
its representation in such a rotating frame, before applying 
the PO procedure, will yield a single condensate mode. We 
consider now a simulation with E = l.lO^g, where in equi- 
librium a single vortex precesses with mean displacement of 
7v ~ 2.0ro from the trap center. After transforming the field 
coeflftcients Cni to a frame (the measurement frame) with fixed 
rotation frequency we construct the covariance matrix in 
the same manner as described in Sec. IV B 1 In Fig. [5ja-b) 
we plot the two most highly occupied modes obtained from 
this covariance matrix, having performed the PO procedure 
in a frame with Qm = 0.2610^;^. We discuss below the rea- 
sons for this choice of rotation frequency. The most highly 



occupied mode is now a vortical state in which the vortex is 
displaced oflT-axis, similarly to the case for the full classical 
field. The second most highly occupied mode is now a peaked 
function which partially 'fills' the density dip associated with 
the vortex in xo- The prominent peak of this mode exhibits 
a distinctive amplitude and phase pattern, which we identify 
as that of the Goldstone mode associated with the breaking 
of rotational symmetry by the oflT-axis vortex state. Such a 
mode arises in scenarios in which the condensate mode 0o 
breaks a symmetry of the system Hamiltonian. In the case of 
broken rotational symmetry |84|, the Goldstone mode has the 
form {uq.vq) ~ (L^0o, -^I^Pq)^ ^nd corresponds to small ro- 
tations of the symmetry-broken condensate mode. As such a 
Goldstone mode is of a fundamentally diflferent nature to the 
'proper' normal modes (quasiparticles) of the system | 85 1, we 
regard this mode as separate from the thermal or nonconden- 
sate fraction of the field |86|. We therefore form the thermal 
density as the total density of all modes other than the con- 
densate mode (xo) and Goldstone mode (xi), 



^th(x) = Y^ni\Xi(x)f. 



(25) 



i>2 



We plot this density in Fig. (Sjc). We note that the thermal 
density exhibits a sharp peak at the vortex location, resulting 
from the erratic motion of the vortex about its mean trajectory 
in the classical-field evolution. Thus a careful application of 
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FIG. 4: (Color online) (a) Mean occupations of covariance matrix 
eigenmodes. (b) Temporal power spectrum of coefficients ai{t) of 
the five most highly occupied modes. Quantities presented in each 
subplot are calculated for the case E - 1.05£g. 



the PO procedure yields a core-filling thermal component, as 
is predicted by the Hartree-Fock-Bogoliubov family of self- 
consistent mean-field theories ISTl . This is in contrast to the 
instantaneous classical field, which, due to its essential na- 
ture as a single-particle wavef unction (i.e., vector in the pro- 
jected single-particle Hilbert space) contains only 'bare' vor- 
tices, with no possibility of core filling. In the microcanonical- 
ergodic approach of the PGPE 131] |35l the two (and higher) 
particle correlations emerge from the field when averaged over 
time. However, an alternative interpretation of classical-field 
trajectories (such as those of |38|) as representative of sin- 
gle experimental trajectories can be made, but, as noted in 
(311, this interpretation must be made with caution. In par- 
ticular it is not correct to interpret the instantaneous state of 
the many-body wavefunction as the Hartree product of the 
classical field ij/ix) (as is the interpretation of the GP wave- 
function in the zero-temperature GP theory [14, 88 1.) Rather, 
the field j/^(x) must be interpreted as an approximation to the 
quantum field operator, from which certain correlations can 



be deduced. Some averaging of this instantaneous field </^(x), 
over an ensemble of trajectories | 48 1 or over time [ 351 must be 
performed to construct the correlations necessary to describe 
the thermal core-filling of a vortex. 

We calculate the time-dependent expansion coefl&cients 
of the modes ;^/(x) in the same manner as described in 
Sec. 



IVBl 



The mean occupations of the modes are plot- 
ted in Fig.[6ja), and reveal that a single mode contains 89.1% 
of the field population. We present in Fig. |5jd) the coher- 
ence functions obtained from the mode coeflftcients. The most 
highly occupied mode is now the only mode with g^^^ signif- 
icantly less than the thermal value g^^^ = 2. We find in fact 
g^Q^ = 1.0002, indicating very high coherence of the conden- 
sate mode over the averaging period. The next most highly oc- 
cupied modes now have values g^^^ > 2.2, greater than that ex- 
pected for thermally or incoherently occupied single-particle 
modes. This is a signature of the strongly collective nature 
of the low-lying excitations of the condensate. The form of 
these excitations is strongly affected by the presence of the 
interacting, macroscopically occupied condensate, which in- 
duces anomalous (pairing) correlations in the field, which in 
turn manifest as anomalous values for the coherence functions 
1891 . We note that the Goldstone mode has g^^^ = 2.90, close 
to the upper bound g^^^ = 3 |89|, as we might expect given 
that a Goldstone mode, with uq = -v^, is the most strongly 
collective excitation possible. The fact that these pair corre- 
lations are manifest in the density matrix formed in the co- 
rotating frame is further evidence that this frame gives the 
correct resolution of the condensate. It is interesting to note 
that the anomalous nature of these correlations is not apparent 
in coordinate- space representations of the field 1351 , or in the 
lab-frame analysis of Sec.|IVBT] A more natural characteri- 



zation of the low-lying excitations would exploit these anoma- 
lous one-body correlations in order to construct an appropriate 
quasiparticle basis, this subject will be discussed elsewhere 
L90|. 

The choice of measurement-frame frequency (Ilm = 
0.26106L)r) optimizes the resolution of the condensate as we 
now discuss. 



IVBl 



we expect 



From the discussion of Sec. 
(neglecting for now the role of the Goldstone mode) that the 
classical field can be resolved into a single condensate mode 
coexisting with a bath of thermal atoms, in some rotating 
frame. In general the rotation frequency of the measurement 
frame diflfers from the equilibrium rotation frequency of the 
condensate, and this diflTerence leads to dephasing of the dis- 
tinct angular-momentum components of the condensate mode, 
and hence to a spurious fragmentation of the condensate into 
a set of temporally coherent angular-momentum eigenmodes. 
We propose as a model for analysis of classical-field data, that 
in the frame in which condensation occurs, the field should 
decompose into a single condensate mode plus thermal mate- 
rial upon performing the PO procedure. This frame depends 
on a single parameter (the measurement-frame rotation fre- 
quency Qm), and the magnitude of the largest eigenvector, 
regarded as a function of the measurement-frame frequency 
[^0 = ^o(^m)] forms the corresponding optimality function. 
Our procedure then is to maximize the optimality function 
to determine the optimal measurement- frame frequency Qc, 



11 




Occupation rii 

FIG. 5: (Color online) Density and phase of (a) the condensate and (b) the Goldstone mode, as determined from the covariance matrix formed 
in the co-rotating frame. The condensate density is the density of mode;^o(x), shown normalized to its occupation hq. (c) Thermal density of 
the field, (d) Second-order coherence functions of the covariance-matrix eigenmodes versus mode occupation. 



which is the frequency of the rotating frame in which the con- 
densate mode is (most) stationary. This provides a sensitive 
measure of the (lab-frame) rotation frequency of the vortical 
condensate mode, and provides a unified basis for the calcu- 
lation of the properties of the condensed and noncondensed 
material in the classical field. In particular, we will regard 
the frame of vortex precession and the frame of condensation 
as interchangeable, though in practice we will use the latter 
exclusively in the remainder of this paper. In Fig. |6jb) we 
present results of this optimization procedure as applied to our 
simulation with E = 1.1 0£g over a 100 eye. period starting 
from t = 9000 eye. We observe a prominent peak in no(^m) 
at Dm = 0.26lO<jJr. We also plot the magnitude of the second 
largest eigenvalue ^i(Qm), and observe that it exhibits a dip at 
this same frequency, indicating the sensitivity of the one-body 
correlations to the frequency Qm- We take the location of the 
peak in ^o(^m) as an estimate for the rotation frequency flc of 
the condensate 1911 . 



3. Temporal decoherence and sample length 

We now consider how the condensate mode determined by 
the above procedure decays as a function of the time over 
which the averaging procedure is performed. As already noted 



in Sec. IV A the phase of the vortex (with respect to any uni- 
form rotation) diflTuses over time. This diflfusion is intimately 
related to the breaking of the rotational symmetry of the sys- 
tem by the presence of the oflT-axis vortex, as we now dis- 
cuss. The condensate orbital formed in the classical field is 
an oflT-axis vortex mode, which is not an eigenstate of angu- 
lar momentum and thus not rotationally invariant. This breaks 
the rotational [S0(2)] symmetry of the Hamiltonian Eq. ^ 
|92|. A well-known consequence of the breaking of a symme- 
try by a dynamical equilibrium in classical mechanics is the 
appearance of a zero-energy normal mode |93|, and this re- 
sult persists in the quantum theory, where broken symmetries 
are similarly accompanied by so-called spurious or Goldstone 
modes |85j. Such excitations are associated with a collec- 
tive motion without restoring force 185] [gH. In the context 
of classical-field theory, we expect that the erratic evolution 
of the field results in random excitation of this motion, corre- 
sponding to fluctuations of the vortex phase. As there is no 
'restoring force' associated with this excitation, over time the 
phase can drift arbitrarily far from any initial value, invalidat- 
ing the interpretation of the equilibrium as a configuration in 
which energy is equipartitioned over normal excitations of the 
symmetry-broken condensate orbital. As a consequence, the 
decomposition of the field into condensate and nonconden- 
sate by the Penrose-Onsager procedure also deteriorates over 



12 




1 2 



Mode 




0.4 0.6 

Oni [units of Ur 



FIG. 6: (Color online) (a) Occupations of the 10 most highly occu- 
pied eigenmodes of the covariance matrix determined in the frame 
with Qm = 0.26 lOtc;^. (b) Dependence of mode populations on the 
rotation frequency of the frame in which the covariance matrix is 
constructed. The inset shows the behavior of uq close to its maxi- 
mum. 



time. This is of course consistent with our expectation that as 
the averaging time becomes large, the statistics we measure 
converge to those of the ergodic density of the microcanoni- 
cal system, which must reflect the rotational symmetry of the 
Hamiltonian |49|. In Fig. [TJa) we plot the values of the two 
largest eigenvalues of the covariance matrix formed by av- 
eraging for varying lengths of time, each average starting at 
t = 9000 cycles. For each averaging time we have optimized 
the condensate fraction over the frequency of frame rotation, 
which causes only small variation in the optimal frequency 
(< 0.003fiojr). We observe that the largest eigenvalue gets 
smaller as the averaging period increases, while the second 
largest eigenvalue gets larger. Over longer periods [inset to 
Fig. [7ja)] the two eigenvalues approach the values obtained 
from the density matrix constructed in the laboratory frame 
(dashed lines), indicating that the orientation of the conden- 
sate mode has difl'used to such an extent that the density ma- 



trix describes a fragmented condensate in all uniformly ro- 
tating frames. In Fig. [TJb) we plot the standard deviation of 
the vortex phases as a function of time, in a frame rotating 
at Q = 0.261066)^. We observe that for the first ^ 75 trap 
cycles the standard deviation remains reasonably steady at a 
small value cre^ ~ 0.05 radians. After this the standard de- 
viation begins to increase dramatically, and this is ultimately 
reflected in the measured eigenvalues [Fig. [TJa)]. We note 
that the standard deviation we have presented is simply that 
of the distribution of vortex phases measured in a single tra- 
jectory, as a function of the length of time over which vortex 
phases are accumulated. To properly characterize the rate of 
rotational difl'usion would require the calculation of the vari- 
ance of the vortex-phase change over an ensemble of similarly 
prepared classical-field trajectories (c.f. [95 1), a numerically 
heavy task which we do not pursue here. We note further that 
the initial period of comparative stability of the vortex phase is 
not a generic feature of the system, in general the excursions 
of the vortex phase relative to its mean precessional motion 
occur unpredictably. We note however that the condensate 
fraction obtained by averaging for only 10 eye. already agrees 
with the values throughout this period of stability to within 
2 - 3%, in agreement with the findings of |35 1. We therefore 
adopt the following approach to determining the condensate 
in the remainder of this paper: we form the covariance matrix 
pij from samples over a period ofT = 10 eye, in a frame 
rotating at frequency We vary this frequency so that the 
largest eigenvalue of the corresponding covariance matrix is 
maximized. We take the frequency of this optimal frame as 
an estimate of the angular velocity of the condensate, and the 
largest eigenvalue as an estimate of the condensate fraction, 
and similarly for subsidiary quantities obtained from the PO 



decomposition (see Sec. |IV Q . We do this for 100 consecu- 
tive periods of 10 cycles, over the period t = 9000 - 10000 
eye. and average the estimates for Dm and /c. In this way we 
sample the short- time dynamics which determine the correla- 
tions of interest, while exploiting the ergodic character of the 
classical-field evolution over longer times. For the simulation 
with E = 1.10£^g, we find from this procedure Qc = 0.2598^^)^ 
and fc = 0.9110, in reasonably close agreement with the val- 
ues obtained from a single averaging over a longer time pe- 
riod. We conclude this discussion by considering the inter- 
pretation of the condensate fraction we measure, in light of 
our abandonment of formal ergodic averaging. We interpret 
the results of the short-time averaging procedure in terms of 
a symmetry-broken representation 1 69 1 of the many-body sys- 
tem. The one-body density matrix we calculate corresponds to 
a particular (mean) orientation of the condensate mode (vor- 
tex phase ^v). Labeling this density matrix by po^, the full 
one-body density matrix is obtained by integrating over all 
possible vortex phases p^^^ = (l/N)f dOy p^^^ |69, 73], with 
a normalization factor. This one-body density matrix does not 
exhibit condensation in the Penrose-Onsager sense, and this 
is true of the physical A^-body system also. Nevertheless, the 
individual one-body density matrix po^ exhibits a condensate, 
and we interpret the correlations it describes as characteriz- 
ing the behavior of a corresponding physical system observed 
to have a particular vortex orientation. With this in mind, we 
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FIG. 7: (Color online) (a) Largest eigenvalues (no and of the 
covariance matrix as a function of averaging time. Inset shows the 
behavior of hq and ni (upper and lower dotted lines, respectively) 
over longer times, where the eigenvalues approach their values ob- 
tained in the laboratory frame (dashed lines), (b) Standard deviation 
of vortex phases measured in a frame rotating at = 0.26 lOo;^. 



will continue to refer to the largest eigenvalue of the density 
matrix po^ constructed from the classical field trajectory as the 
condensate fraction. 



C. Rotational properties of the field 

The prescription we have developed for decomposing the 
classical field into condensed and noncondensed parts allows 
us to extract quantities which characterize the rotation of the 
two components individually, by calculating appropriate ex- 
pectation values. We recall that the extension of a single-body 
operator / to a second-quantized operator J = f drifr'^Jifr on 
the many-particle Fock space 1961 has expectation value in 
the many -body state (/) = Tr{p^^^/}. We thus define an ex- 
pectation value of a single-body operator in the classical field 



</)e ^ Tr{p/}, 



(26) 



where p is the covariance matrix introduced in Sec. IV B 



The decomposition into a condensed mode and noncondensed 
modes introduced in Sec. |IVB^ [97] allows us to write (us- 
ing Dirac notation for vectors in the projected single-particle 
space) 

P = ^oko)C\rol + ^ nk\Xk)<^k\ = Po + Pth. (27) 

k>0 

We thus define averages in the condensate and noncondensate 
by 

</)o = Tr{po/} = no Yj(^\Xo}<Xo\ J\m}, (28) 



and 



(Ah = Tr{pth/} = Yj^k Yj<'^\^k}<Xk\J\m\ (29) 



k>0 m 



respectively. It is clear from the form of po and pth that they 
are (proportional to) mutually orthogonal projectors, and so 
we have the additivity property (/)o + (/)th = {J)c' 



1. Angular momentum of the condensate 



As noted in Sec. IV B 2 



in the preces sing- vortex sce- 
nario the condensate determined by our method is a non- 
axisymmetric state with an off-center vortex. Analogous con- 
densate orbitals appear in the zero-temperature GP theory as 
intermediates between the rotationally invariant £ = and 
£ = I modes, with fractional angular-momentum (per atom) 
expectation values < ^ < 1 1 16, 98|, and are mechanically 
unstable (see the discussion in |99 |). 

Applying the above averaging procedure to our simulation 
with E = l.lO^g we obtain = (L^^/^o = 0.63^. We re- 
call that the condensate fraction obtained for this state was 
/c = 0.911. Turning to the noncondensate we find Lth/Mh = 
4.80li. While at this energy, only 9% of the atoms have been 
excited out of the condensate, in order to maintain rotational 
equilibrium of the system, the thermal atoms carry nearly half 
of the angular momentum of the field. This is a consequence 
of the suppressed moment of inertia of the condensate: in or- 
der to have the same angular velocity, the thermal cloud must 
possess much more angular momentum per particle than the 
condensate ifTOOlfToTll . 



2. Angular velocity of the thermal cloud 

Now we consider the angular velocity of the thermal cloud. 
In 1 38 1 we gave an approximate method of separating the ro- 
tational properties of the condensate and the thermal cloud, 
based entirely on spatial location. In the current paper, where 
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we have developed a rigorous procedure for separating the 
condensate and thermal components, we can give a more ac- 
curate analysis of the rotation and inertia of both of these two 
components. In particular our approach includes the contri- 
bution of thermal fluctuations that traverse the central con- 
densed region, as well as those which fill the vortex core, to 
the total thermal component of the field. We thus calculate 
the averages (L^)th and (©c)th, and estimate the cloud rotation 
frequency by 



'z/th 



<©c)th 



(30) 



i.e., we assume that the thermal component's moment of in- 
ertia is equal to its classical value. Doing so we find Qth = 
0262(jJr, in fair agreement with the precession frequency of 
the vortex. 



E = 1.04£'g to nearly zero by £" = 1.16£'g. This is a sig- 
nature of the anomalous rotational properties of the off'-axis 
vortex state |98 |: as the energy is increased the condensate 
loses angular momentum to the cloud and the rotation rate of 
both components increases. 

At higher energies the condensate is vortex-free and thus ir- 
rotational, and exists in thermal and diff'usive equilibrium with 
the thermal cloud which now contains all the angular momen- 
tum of the field. The small residual angular momentum of 
these high-energy condensates shown in Fig.[9ja) results from 
surface excitations which are not diff'erentiated from the con- 
densate by the short- time averaging. We note however that the 
condensate fraction is essentially continuous across this tran- 
sition from a symmetry-broken condensate orbital to a non- 
symmetry-broken (i.e., vortex-free) one, supporting the defi- 



nition introduced in Sec. IV B 3 as an appropriate measure of 
condensation in the symmetry -broken regime. 



V. DEPENDENCE OF SYSTEM STATE ON INTERNAL 
ENERGY 



C. Condensate and cloud rotation rates 



We now consider the eff'ect of varying the energy of the 
classical field on its equilibrium properties. We expect quite 
generically for Hamiltonian classical-field simulations such as 
those performed here that an increase in the classical-field en- 
ergy will result in an increase in the field temperature and a 
suppression of the condensate fraction |[35l . In this scenario 
with finite conserved angular momentum, we also expect the 
rotational properties of the system to depend strongly on the 
field energy. 



A. Density 

The most obvious and generic consequence of increasing 
the energy of the microcanonical system is an increase in the 
entropy of the system. The qualitative eff'ects of this increase 
can be readily observed in position-space representations of 
the classical field density. In Fig.[8ja-f) we plot representative 
densities of the classical field at equilibrium, for various val- 
ues of the field energy. These densities are chosen at times in 
the range t e [9000, 10000] eye. such that the vortex is dis- 
placed along the x axis, for ease of comparison. We clearly 
observe from these images an increase in the surface excita- 
tions of the condensate and the development of a turbulent 
outer cloud as the energy is increased. The images also re- 
veal an increase in the vortex precession radius as the energy 
is increased. 



B. Condensate fraction and angular momentum 



In Fig. |9jb) we plot the condensate rotation frequency, as 
determined by the frame frequency of optimal condensate res- 



From the procedure introduced in Sec. IV B 3 we find the 
condensate fraction [Fig.[9ja)] drops gradually from /c = 0.96 
to fc = 0.83 over the energy range E e [1.04, 1.19]£^g. 
However, the angular momentum (per particle) of the con- 
densate drops dramatically, falling from (L^)o/^o = 0.86^ at 



olution (Sec. |IVB2| ). We find that the rotation frequency in- 
creases gradually from = 0.234^^)^ at £ = 1.04£^g, up to 
Qc = 0.340^;^ at E = l.lSEg. For comparison the frequency 
of the anomalous mode of the ground vortex state in the Bo- 
goliubov approximation, kanoml/^ = 0.20346l)^, is plotted as 



a dotted line in Fig.[9jb). We note that the case E = 1.15£^g 
is very close to the transition from a vortical condensate to 
a vortex-free one; in ~ 10% of the 10 eye. sampling peri- 
ods the condensate is vortex-free, resulting in maximal val- 
ues of ^o(^m) occurring around Qm = l^r |102|. We dis- 
card these values in calculating the mean and standard devia- 
tion presented for this simulation in Fig.[9jb). From Fig.[9ja) 
we observe that the condensate modes with E > 1.1 6£g are 
essentially irrotational, and indeed the condensates in these 
simulations appear to be largely vortex-free, the majority of 
samples yielding Qc ~ l^r, with a small minority produc- 
ing Dc ~ OAoJr, resulting from the transit of short-wavelength 
surface excitations (so-called ghost vortices) which produce 
a slightly greater calculated condensate fraction in a frame 
which matches their motion, and we omit these points from 
Fig.igb). 

Turning our attention to the cloud rotation frequency 
[crosses in Fig.jg^b) (red online)], determined by the prescrip- 
tion of Sec. |IVC2l we find these results in qualitative agree- 
ment with the measured condensate rotation frequency over 
the range E e [1.04, 1.15]£'g. The bars accompanying the 
markers for condensate and cloud rotation rates simply rep- 
resent the statistical dispersion in the measured quantities as 



discussed in Sec. |IVB 3] and not quantitative measures of the 
errors in these measurements, but it is clear that the measured 
cloud rotation rate does not match the condensate rotation rate 
closely. It is possible that this is associated with pair correla- 
tions in the noncondensate, though if anything we expect the 
presence of pair correlations to suppress the moment of iner- 
tia relative to its classical value (c.f. 11 03 J ). We note finally 
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that the cloud frequency reaches a maximum of Qth ~ 0.33^^^ 
at E = 1.16£^g before decreasing as the energy is increased 
further. This reduction of Qth with increasing energy is easily 
understood: once the condensate becomes vortex-free (irrota- 
tional), as further atoms are depleted from the condensate the 
thermal cloud population increases while its angular momen- 
tum remains fixed, implying a decrease in angular velocity. 



D. Temperature and chemical potential 

In order to estimate the temperature and chemical poten- 
tial of the field, we fit the distribution of thermal atoms in the 
classical field with an approximate semiclassical distribution 
(381. The principle of the method is to form a fitting func- 
tion corresponding to the semiclassical distribution of atoms 
in the combined potential formed from the centrifugally di- 
lated trapping potential and the mean-field repulsion of the 
(time-averaged) classical field. In contrast to |38|, the rota- 
tion rate of the thermal cloud, which defines the centrifugal 
dilation of the trapping potential experienced by the thermal 
atoms, is not assumed to be the same as that of the frame in 
which the projector is applied. We are therefore obliged to 
generalize the fitting function of 1 38 1 to allow for a differential 
rotation between the thermal cloud and projector. The details 
of this generalization are included in Appendix |B] where we 
derive the form of the fitting function 

n(nji,T) = -^[/i(r;//,Dth,np) + /2(r;yu,Dth,^^p)]. (31) 



In order to avoid including significant condensate density in 
our fitting, we are obliged to restrict our fit to a domain 
re [r_, rtp], where r_ is the location of the minimum of the 
eff'ective potential, which we expect to approximately mark 

the condensate boundary, and rtp = ^2ERlm{ojl - Dp) is 
the semiclassical turning point of the condensate band in the 
frame of the projector |[38]| . In practice the low tempera- 
tures of the classical-field configurations we consider make 
performing such a fit difficult, in particular the fit fails close 
to r_, in the region of strongly collective excitations of the 
condensate. We therefore perform our fits only to the wing 
of the classical-field density; for definiteness we take r e 
[r_ -h Iro, rtp]. We take the values Qth calculated in Sec. VC|as 



estimates for the cloud rotation rate in Eq. (31). Although 



these values may be somewhat inaccurate as discussed in 



Sec. VC we find that using the vortex precession frequen- 
cies (where available) makes little diff'erence to the obtained 
values of jj, and T, and conclude that any discrepancy in Qth 
is inconsequential at the level of accuracy of the semiclassi- 
cal estimates we seek. The results of this fitting procedure 
are presented in Fig. [9jc). We observe that the temperature 
estimate RbT increases approximately linearly with the en- 
ergy E, and is of the same order of magnitude as the chemical 
potential yu. That the chemical-potential estimates are gen- 
erally higher than the chemical potential of the ground state 
yUg = 10 35fiojr, whereas we expect the chemical potential 
of the strictly Hamiltonian system to decrease with increas- 
ing energy (see 1 104|) suggests that this quantity is probably 
overestimated by the fitting, and some error in the tempera- 
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FIG. 9: (Color online) Dependence of equilibrium parameters on field energy, (a) Condensate fraction and condensate angular momentum, 
(b) Angular velocities of the condensate and thermal cloud, (c) Temperature and chemical potential, (d) Vortex precession radius, classical 
moment of inertia of the thermal cloud, and nonclassical moment of inertia of the condensate mode. Errorbars indicate the statistical spread in 
values calculated (see main text). 



ture estimate is also expected. Nevertheless the results suggest 
that our simulations probe the temperature regime > jdin 
which the thermal friction on the vortex is expected to be ap- 
preciable 1 27 1, consistent with our observations. 



imposed by the thermal cloud, and its moment of inertia is 
thus zero. 



E. Moment of inertia and vortex displacement 

In Fig. [9jd) we plot the classical moment of inertia of the 
noncondensate (crosses), <0c)th (see Sec. |IVC 2| ). We find that 
this quantity increases steadily with the field energy as atoms 
are excited out of the condensate into the thermal cloud. We 
also plot the (quantum) moment of inertia of the condensate 
orbital (plusses), which we define as ©o = (^z)o/^c, where 
the angular momentum (L^)o and rotation frequency Qc of 
the condensate are those discussed in the previous sections. 
To give an indication of the variation of estimates in this de- 
rived quantity, we include the standard deviation estimated by 
adding the standard deviations of (L^)o and Dc in quadrature 



I1Q5I . We note that the moment of inertia of the condensate 
decays steadily with increasing E, following (L^)o, and by 
E = 1.15£^g is essentially zero. Above this energy the conden- 
sate's angular momentum is basically zero (see the remarks in 



Sec. VB ), i.e. the condensate fails to respond to the rotation 



On this plot we also include the displacement of the vortex 
from the trap center (circles). For consistency with the other 
measurements we calculate this by tracking vortex trajecto- 
ries over 100 successive intervals of 10 eye. (starting from 
t = 9000 eye.) and averaging the measured radii. Bars on 
the marker again indicate standard deviations of the estimates 
over the 100 periods. The precession radius increases steadily 
as the energy is increased. By £" = 1.13£g, the precessing 
vortex is close to the edge of the condensate, and so the sim- 
ple vortex tracking algorithm begins to have difl&culty tracking 
this single vortex in proximity to the proliferation of phase 
defects undergoing constant creation and annihilation in the 
condensate periphery |38|. We therefore do not attempt to in- 
clude data for these higher energies on the plot. It is clear 
however that the vortex displacement increases as the conden- 
sate angular momentum decreases, as is the case in the zero- 
temperature GP theory |[T6l[T9HT06ll . 
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VI. SUMMARY AND CONCLUSIONS 

We have carried out classical-field simulations of a precess- 
ing vortex in a finite-temperature Bose-Einstein condensate. 
Employing a microcanonical ergodic-evolution approach, we 
find that randomized states of the classical field subject to an 
additional angular-momentum constraint settle to configura- 
tions each containing a single vortex precessing under the in- 
fluence of its own induced velocity field, at thermal and rota- 
tional equilibrium with the thermal cloud, so that the mutual 
friction force between the vortex and the cloud vanishes. 

As the classical field condenses into a rotationally 
symmetry-broken state, a naive application of the Penrose- 
Onsager test of one-body coherence in an inertial frame yields 
a condensate fragmented over a set of coherent angular- 
momentum eigenmodes. Transforming to a rotating frame 
such that the fragmentation is eliminated we find that the clas- 
sical field describes a condensate mode with an off'-axis vortex 
coexisting with a core-filling thermal component, and we also 
observe the Goldstone mode associated with the rotational 
symmetry breaking. 

Due to the ergodic nature of the classical-field evolution, 
the rotational phase of the vortex diff'uses over time. Thus care 
must be taken to separate the short-time fluctuation dynamics 
which define the coherence of the field from the dynamics on 
much longer time scales, over which ergodic migration of the 
field configuration occurs. We found that extracting the con- 
densate by averaging over short times and averaging over suc- 
cessive estimates so obtained allowed us to access the short- 
time dynamics of interest while exploiting the ergodic nature 
of the system on longer time scales to better estimate system 
observables. 

We showed that the angular momentum of the condensed 
mode decreases dramatically as the energy (and thus tempera- 
ture) are increased at fixed total angular momentum, and cor- 
respondingly the vortex precessional radius and frequency in- 
crease until the vortex leaves the condensate, and the conden- 
sate becomes irrotational. We proposed a measure of the cloud 
angular velocity based on the decomposition of the field into 
condensed and noncondensed components and found qualita- 
tive agreement with the vortex precession frequency. 

Possible extensions of the work presented here are to 
nonequilibrium scenarios of vortex decay and to multiple vor- 
tex configurations. We also expect the approaches developed 
here to be of more general applicability in extracting dynam- 
ical characteristics of condensates in ergodic classical field 
studies. 
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FIG. 10: (Color online) Dependence of equilibrium system param- 
eters on the frame in which the projection is effected, (a) Angular 
velocity of the condensate, (b) condensate fraction, (c) vortex pre- 
cession radius, (d) temperature and chemical potential. 



APPENDIX A: DEPENDENCE ON FRAME OF PROJECTOR 



We now consider the effect of the angular velocity flp 



of 



the frame in which the projector is defined on the equilib- 
rium configurations of the classical-field trajectories. In ir- 
rotational scenarios in which the projector is applied in a lab- 
oratory frame, the cutoff energy has a significant effect on the 
equilibrium attained for given conserved first integrals (en- 
ergy and normalization). In the present case, the cutoff is de- 
fined by two parameters: the angular velocity of the frame in 
which the cutoff is effected, and the (rotating-frame) energy 
at which the cut is made. As noted in Sec. [Ill B 2| the (ther- 
modynamic) angular velocity of the classical field at equilib- 
rium is in general not equal to that of the projector. Here we 
consider the effect of the choice of projector angular veloc- 
ity on the classical-field equilibrium. We vary the projector 
frequency over the range - 0.35 cl)^, and keep the (rotating- 
frame) cutoff fixed at Er = 30fiajr. There are several conse- 
quences of varying the angular velocity Qp. Firstly the mul- 
tiplicity At (i.e., the number of canonical coordinate pairs in 
the Hamiltonian system) depends strongly on Qp, increasing 
with increasing Moreover the set of modes included varies 
with Dp, with the condensate band becoming increasingly bi- 
ased towards single-particle modes with m > as Dp is in- 
creased 1 107 1 . In Fig. [TO|a-d) we present the results of sim- 
ulations with E[i//] = l.lO^g, Er = 30nojr, L[i/^] = NqK and 
Qp G [0, 0.35]^;^. We find that all these simulations exhibit 
qualitatively similar equilibrium states, each containing a sin- 
gle precessing vortex, as for the simulation with E = 1.10£^g 



discussed in the main text. In Fig. [TO[a) we plot the angular 
rotation frequency of the condensate mode, as determined by 



the procedure of Sec. IV B 3 The rotation frequency exhibits 
a clear upward trend as the frame frequency is increased, nev- 
ertheless the dependence of the rotation rate on frame fre- 
quency is weak, varying by only ~ 0. 03 ojr over the range 



of Qp considered. In Fig. 10 ^b) we plot the corresponding 
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condensate fractions. Again the dependence on the frame fre- 
quency is weak; any trend in the condensate fraction is of the 
same magnitude as the statistical dispersion of the measured 
values. In Fig.fTOtc) we plot the vortex precession radii. We 



see that the radius increases steadily, approximately doubling 
between Ilp = and Qp = 0.35^;^. This increase in preces- 
sion radius corresponds to a loss of angular momentum (per 
particle) from the condensate, and is consistent with the ob- 



served increase in condensate angular velocity. In Fig. \lU[d) 
we plot the dependence of the chemical potential (crosses) and 
temperature (circles) on Dp, calculated assuming the conden- 
sate angular frequency Dc as an estimate for the cloud rota- 
tion rate. We observe that the temperature decreases as the 
frame frequency is increased, consistent with the increase in 
the number of modes in the condensate band (c.f. |38 1). We 
conclude then that increasing the frame angular velocity Qp 
results in a lowering of the equilibrium temperature, and bi- 
ases the angular-momentum balance in the system towards the 
thermal cloud, leading to a higher equilibrium angular veloc- 
ity. We emphasize that this dependence is a natural conse- 
quence of the projector definition: diff'erent values of Dp de- 
fine diff'erent Hamiltonian systems. Ultimately, this depen- 
dence would be removed by the inclusion of above-cutoff' ef- 
fects 13111331 [361. 



APPENDIX B: SEMICLASSICAL FITTING FUNCTION 

In order to estimate the thermodynamic parameters of the 
field, we assume an approximate semiclassical description for 
the thermal atoms. The description consists of two elements: 
a classical-field distribution in the phase-space coordinates (x 
and p) 1 38 1, and the specification of an appropriate energy 
cutoff'. In order to construct this distribution, we must take 
into account two generally disparate rotations: the equilibrium 
rotation of the thermal atoms (angular velocity Ilth) and the 
rotation of the projector defining the cutoff' (angular velocity 



Qp). The classical-field distribution takes the general form 



^c(x, P) = 



kBT 



6(X, p) - yU 



(Bl) 



where the semiclassical energy of the thermal atoms (rotating 
at Qth) is 



6(X,P)=^ + 

2m 



2 2 



- Q.ih{xpy - ypx) + 2An(r), (B2) 



(c.f. 11081 ) and where n(r) is the total density of the atomic 
field, which we will take to be circularly symmetric. The ap- 
plication of the projector P defined in a rotating frame en- 
forces a cutoff' in the energy, which delimits the accessible 



region of phase space occupied according to Eq. (Bl ). This 
region is defined by 



£1 
2m 



2 2 

mcotr 



■ Qp(xpy - ypjc) + 2A.n{r) < Er. 



(B3) 



2 -^p\-^Fy 

The density of thermal atoms at any spatial position x is 
obtained by integrating Eq. (Bl) over the momentum, i.e. 
^(x) = ^ JpFc(x, p), where the domain D in momentum 
space is determined by Eq. (B3), and is itself a function of 
the position-space coordinate x. To perform this integration 
we make the transformation to the kinematic momentum mea- 
sured in the frame of the cloud P = p - mOth x x l l08[ |lQ9ll . 
In this representation the integrand becomes circularly sym- 
metric in P, i.e., we now have F^(x,P) = k^T l[e(x,V) - yu], 
where 



e(x,P) = — + -K - ^IV + ^Mr). 



(B4) 



In these momentum coordinates the domain D is circular, 
however its center is in general displaced from the origin 
P = 0. After some simple algebra we find that the domain 
of integration is 



l)(x) : —{[P, + m(Qp - Qth)j]' + [Py ■ 



■m(^^-QM^\<ER- 



-micOj. 



■ - 2An{r). 



(B5) 



As the integrand is circularly symmetric in P, the result 
of the integral depends only on the radius r of the chosen 
point in position space, as azimuthal rotations of the coor- 
dinate X simply induce rotations of D(x) about the origin 
of momentum space. We therefore consider without loss of 
generality {x,y) = (0, -r). The domain D thus becomes 
a circle centered on Px = m(Qp - Qth)^ = a, of radius 



l2m[ER - (m/2)(cof - Q^)r2 - 2An(r)]. Depending on 

the value of r, we may have a < h or b > a. We adopt polar 
coordinates (P, 0), and find for the angle A0 subtended by D 



at a given radius P 



A(P(P) = 



2 cos 
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-1 ^p^ 





2aP 



P < a - b (when a > b) 
P < b - a (when a < b) 
) \a-b\<P<a-\-b 
P > a + b 



(B6) 
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Integrating in polar coordinates we obtain finally 
27imkBT PdPA(p(r) 



n(r) 



p PdPA 

Jo £7^^^ 



2An(r) - [i 



4b 



^ [/i(r; Dth, ^^p) + /2(r; ^th, ^^p)], (B7) 

where 

/i(r) = &{b-a)\n 
and 



Ka;2-Q2^)r2 + 2Mr)-A/ 



(B8) 



/2(r) 



+^ PdP cos-^ ( 



1 ( ^ + 'f[ajj-Q.lHnp-Q.,hfy-ER 



(B9) 



and we have introduced the de Broglie wavelength A^b = 
^J27:fi^/mkBT |110|. The integral I2 must be evaluated nu- 
merically. In the limit ^th ^p the integral I2 vanishes, and 
we regain the form 



n(r) 



employed in 



In 



Er-ji 



{c4-al)r^ + 2An{r)-ii 
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